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1.  Introduction 


Improved  accuracy  of  fire  and  efficient  design  of  anti¬ 
recoil  devices  cannot  be  achieved  without  a  detailed  knowledge  of 
the  flow  field  in  the  muzzle  region  of  a  gun,  in  the  short  inter¬ 
val  of  time  (of  the  order  of  a  few  milliseconds,  at  most)  between 
the  first  appearance  of  the  projectile  and  its  exit  through  the 
precursor  shock  wave.  Experimental  techniques  have  evolved  to  a 
high  degree  of  sophistication  [1,19,20,21]  so  that  a  complete 
description  of  one  firing  can  be  obtained  by  a  combination  of 
measurements,  photographs  and  data  processing.  Innovative 
design,  however,  requires  a  similar  detailed  analysis  of  a  large 
number  of  cases,  using  different  models;  the  task  is  too  time- 
consuming  and  expensive  to  be  accomplished  experimentally.  Nu¬ 
merical  analysis  seems  to  be  a  possible  alternative  to  experi¬ 
ments,  although  results  prior  to  the  present  work  have  not  been 
too  encouraging. 

The  numerical  techniques  used  so  far  have  been  largely 
inadequate  to  the  task.  Certain  requirements  must  be  satisfied, 
indeed,  viz.: 

1)  The  geometry  of  the  muzzle  and  of  the  projectile  must 
be  realistic, 

2)  The  computational  program  must  be  able  to  handle  more 
complicated  geometries  than  just  a  muzzle  and  a  projectile  (for 
example,  it  should  include  deflectors  of  arbitrary  forms), 

3)  The  motion  of  the  project’ le  must  be  realistic, 

4)  The  boundary  conditions  at  all  rigid  walls  must  be 
properly  handled, 

5)  No  artificial  outer  boundary,  capable  of  introducing 
systematic  numerical  errors,  should  exist  in  the  program, 

6)  Shocks  should  be  handled  as  discontinuities  satisfy¬ 
ing  the  Rankine-Hugoniot  conditions, 

7)  Prediction  of  imbedded  shock  formation  and  treatment 
of  shock  interactions  should  be  included  in  the  code, 

8)  Two  different  gases,  having  different  molecular 
weights  and  different  ratios  of  specific  heats,  should  be  con¬ 
sidered  (although  their  thermodynamical  behavior  could  be  assumed 


Introduction 


to  be  that  of  a  perfect  gas), 

9)  Consequently,  contact  discontinuities  should  be  ex¬ 
plicitly  handled  by  the  code, 

10)  Realistic  initial  conditions  should  be  used. 

Practically,  none  of  the  requirements  listed  above  i3  sa¬ 
tisfied  by  any  of  the  existing  codes  [2-9].  By  and  large,  most 
of  such  codes  are  inspired  by  the  early  work  of  the  Los  Alamos 
group  [10-12]  which  has  undisputed  pioneering  merits  but,  in  the 
light  of  modern  developments,  appears  to  be  representative  of  a 
'brute  force'  approach,  typified  by  the  use  of  Cartesian,  evenly 
spaced  grids,  rudimentary  boundaries  running  along  grid  lines, 
simplistic  treatment  of  boundary  conditions,  smearing  of  discon¬ 
tinuities  and  low  order  of  accuracy  of  the  integration  schemes. 

A  cursory  glance  at  the  computational  grids  used  in  the 
Reports  cited  above  shows  not  only  their  inadequacy  to  provide 
appropriate  resolution  where  needed  but  the  unrealistic  geometry 
of  rigid  walls  as  they  result  from  the  choice  of  a  basic  Carte¬ 
sian  mesh  (in  other  words,  the  mesh  determines  the  wall  geometry 
whereas  the  opposite  should  occur  in  a  realistic  computation). 
In  handling  rigid  boundaries,  large  use  of  reflective  cells  is 
made  (a  by-product  of  Cartesian  grids  and  straight  boundaries); 
this  can  be  a  major  source  of  errors  [13].  The  outer  boundaries 
of  the  computational  region  are  generally  mishandled  as  well. 
Only  recently  ha3  the  majority  of  numerical  analysts  accepted  the 
fact  that  disturbances  produced  inside  the  computational  region 
must  be  allowed  to  propagate  outwards  without  reflecting  on  any 
artificial  outer  boundary.  How  to  achieve  such  a  goal  in  prac¬ 
tice,  however,  is  not  quite  clear  yet.  SAMS  (the  most  advanced 
code  to  date  and  the  only  code  which  reflects  a  real  concern  for 
the  treatment  of  boundary  conditions)  sets  all  gradients  to  zero 
at  the  outer  boundaries.  Such  a  procedure  is  inconsistent  with 
the  outwards  physical  propagation  of  disturbances  through  the 
boundaries.  Thi3  mistake  was  noticed  by  Zoltani  [14]  who  tried 
using  semi-empirical  mass  sinks.  We  believe  that  this  approach 
i3  not  correct  and  might  be  responsible  for  some  of  the  catas¬ 
trophic  oscillations  which  the  modified  SAMS  code  generates  and 
which,  according  to  Zoltani,  are  not  smoothed  by  standard  artifi¬ 
cial  viscosity  devices.  Finally,  no  fitting  of  shocks  and  con¬ 
tact  discontinuities  has  been  attempted  so  far;  to  avoid  numeri- 
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cal  instabilities,  some  forms  of  damping  have  to  be  provided, 
which  deface  the  pattern  of  shocks  and  contact  discontinuities, 
typical  of  the  nature  of  the  muzzle  flow. 

The  analysis  presented  in  this  paper  has  been  developed 
on  a  philosophy  totally  different  from  any  previous  approach. 
Instead  of  trying  to  solve  the  problem  on  the  whole  field  of  in¬ 
terest  considered  as  a  single  computational  region,  and  using  a 
single,  general-purpose  integration  technique  for  the  equations 
of  fluid  motion  over  the  entire  area,  we  subdivide  the  region  of 
interest  into  a  number  of  subregions,  the  inner  and  outer  boun¬ 
daries  of  which  are  defined  on  physical  grounds.  Each  region  can 
be  studied  in  detail;  accuracy  is  easier  to  achieve  and  informa¬ 
tion  i3  not  lost.  In  addition,  the  analysis  of  each  region  is 
based  on  the  most  sophisticated  and  efficient  numerical  tech¬ 
niques  currently  available,  namely,  the  use  of  conformal  mappings 
to  define  computational  grids,  the  use  of  the  X-scheme,  the  in¬ 
clusion  of  viscous  effects,  and  the  fitting  of  shocks  and  contact 
discontinuities. 

In  the  present  paper,  we  describe  the  calculation  from 
the  instant  of  firing  through  the  instant  at  which  the  bullet  to¬ 
tally  emerges  from  the  barrel.  The  analysis  is  divided  into 
three  parts,  originally  handled  by  separate  codes: 

1)  One  dimensional  flow  inside  the  barrel,  in  front  of 
the  bullet  (code  name,  NS14) 

2)  Axisymmetric  flow  around  the  muzzle,  with  the  bullet 
still  inside  the  barrel  (code  name,  NS9),  and 

3)  Axisymmetric  flow  around  the  muzzle  during  the  exit 
of  the  bullet  (code  name,  NS18). 

The  three  codes  have  been  combined  into  a  single  one, 
named  NS21 ,  which  provides  an  uninterrupted  description  of  the 
precursor  phase  of  evolution,  based  solely  on  geometrical  and 
internal  ballistic  data. 
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2.  One-dimensional  flow  inside  the  barrel 


In  this  calculation,  the  bullet  is  assumed  to  be  flat¬ 
faced.  This  is  not  a  major  restriction,  due  to  the  scale  of  the 
ogive  with  respect  to  the  length  of  the  barrel,  and  to  the  fact 
that  the  gas  near  the  bullet  front  tends  to  move  a3  a  rigid  body. 
A  law  of  motion  of  the  bullet  as  a  consequence  of  the  firing  of 
the  charge  is  assumed  from  empirical  data.  Generally,  the  bullet 
moves  with  no  initial  acceleration;  consequently,  the  first  per¬ 
turbation  is  just  a  characteristic  moving  into  the  barrel  at  the 
speed  of  sound  of  the  gas  at  rest.  Shortly  thereafter,  an  imbed¬ 
ded  shock  forms  in  the  perturbed  region  and  rapidly  overtakes  the 
perturbation  front.  By  the  time  the  imbedded  shock,  which  Is  now 
identified  with  the  perturbation  front  and  thus  became  the  pre¬ 
cursor  shock,  reaches  the  gun'3  muzzle,  the  flow  between  the 
shock  and  the  bullet  is  almost  uniform. 

The  calculation  is  performed,  at  every  step  in  time, 
between  the  perturbation  front  and  the  face  of  the  bullet.  Times 
are  counted  from  the  instant  the  bullet  starts  moving.  The  cal¬ 
culation  begins  at  a  time,  tQ,  which  i3  very  small.  The  pertur¬ 
bation  front  is  at  a  distance  from  the  initial  location  of  the 
bullet  face  equal  to  aQt0«  where  aQ  is  the  speed  of  sound  of  the 
gas  at  rest.  A  number  of  nodes  is  chosen  to  define  the  flow;  in¬ 
itially,  only  two  intervals  are  considered  between  bullet  and 
perturbation  front ;  the  number  of  intervals  is  doubled  as  the 
distance  between  bullet  and  perturbation  front  exceeds  a 
prescribed  multiple  of  the  initial  distance;  the  doubling  is  re¬ 
peated  over  and  over  again,  u3ing  the  same  criterion,  until  a 
maximum  of  6U  intervals  is  reached.  If  the  distance  decreases, 
provision  is  made  for  halving  the  number  of  intervals. 

The  length  of  the  computational  region  is  normalized; 
that  is,  in  addition  to  the  physical  space  coordinate,  z,  and  to 
the  physical  time,  t,  we  use  the  computational  variables,  X  and 
T,  related  to  z  and  t  by 
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X  =  (z-b)/(c-b)  ,  T  =  t  (1) 

where  b  and  c  are  the  z -coordinates  of  the  bullet  face  and  of  the 
perturbation  front,  respectively. 

The  flow  in  the  barrel  is  assumed  to  be  inviscid.  Let  p, 
p,  S  and  9  be  the  thermodynamical  parameters  density,  pressure, 
entropy  and  temperature,  respectively,  u  the  velocity  and  t  the 
time.  Pressure  and  density  of  the  gas  at  rest  are  assumed  as  un¬ 
ity;  the  speed  of  sound  of  the  gas  at  rest,  divided  by  the  square 
root  of  y,  is  assumed  as  the  unity  of  speed.  The  inner  radius  of 
the  barrel  is  assumed  as  the  unity  of  length.  With  a  suitable 
definition  of  the  unity  of  time,  as  the  ratio  of  the  unity  of 
length  to  the  unity  of  speed,  and  letting 

.  P  =  In  p  ,  9  =  p/p  (2) 

S  =  y  In  9  -  (y-1)  P  (3) 

the  equations  of  motion  in  the  (X,T)  frame  are: 

PT  +  APX  -►  yXzux  =  0 

uT  +  Aux  ♦  6  Px  =  0  (4) 

S.j.  ♦  ASX  =  0 

where 

A  =  uXz  ♦  Xt  (5) 

Xz  =  1 / ( c-b )  ,  Xt  =  X  [(X-l)bt  -  Xct]  (6) 

The  equations  of  motion  are  integrated  at  all  nodal 
points  (including  the  node  on  the  bullet  and  the  one  at  the  per¬ 
turbation  front)  using  the  i-scheme  [15].  Special  treatments  are 
given  to  the  boundary  points  and  the  points  next  to  the  boun¬ 
daries,  to  avoid  using  information  from  outside  the  boundaries. 
Note  that  at  the  bullet  point  we  need  to  compute  the  pressure  on¬ 
ly,  since  u  coincides  with  the  prescribed  b  and  S  is  always 
equal  to  its  initial  value,  assumed  as  zero.  The  pressure  is. 
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thus,  obtained  by  using  the  equation  for  left  running  charac¬ 
teristics  : 


=  -SPX1 


a  (btt" 


Vxi5 


A  -  X  a 
z 


(7) 


At  the  end  of  the  corrector  level,  at  each  time  step,  pressure 
and  velocity  are  corrected  at  the  perturbation  front,  assuming 
that  it  is  a  3hock  (obviously,  in  the  first  phase  of  the  motion, 
such  a  shock  degenerates  into  a  characteristic).  Shocks  are  com¬ 
puted  using  the  post-correction  technique  [16].  In  the  case  of  a 
shock  moving  into  a  gas  at  rest  with  a  speed  W  =  c  ,  we  have 


3P/3W  =  2W/(W  -6) 


3u  .  /3W  =  -  S(6«.1)  ♦  y/[(«*1)W2] 
rel 

i  =  (y-1)/2 


(8) 


P*  =  ln(W2-5)  -  ln(«+1)  ,  u*  s  -  <S/(4*-1)W  -  r/[«*DW]  (9) 

rel 

AW  =  -(a(P*-P{E))+Y(u*  ,+W-u(E))]/[a  3P/3W  +  y ( 3u  , /3W+1 ) ] ( 10) 

rel  rel 

(E )  ( E ) 

where  P  ,  u  are  the  values  obtained  by  integrating  the  Euler 
equations.  After  updating  W,  the  values  behind  the  shock  are 
recomputed  from  the  Rankine-Hugoniot  conditions.  The  location  of 
the  shock,  c,  is  updated  using 

Ac  =  (ct  +  ^cttAt)At  =  (W  +  ^AW)At  (11) 


At  every  computational  step,  the  value  and  location  of 
the  maximum  pressure  gradient  is  monitored;  when  it  exceeds  a 
given  tolerance  (say,  3P/3z>2),  an  imbedded  shock  is  fitted  at 
the  center  of  the  computational  mesh  where  the  maximum  gradient 
occurs.  The  interpolated  local  pressure  is  assumed  to  be  the 
pressure  in  front  of  the  shock.  The  pressure  at  the  next  node  is 
assumed  as  the  pressure  behind  the  shock;  and  the  shock  Mach 
number,  initial  velocity  and  other  pertinent  values  are  calculat¬ 
ed  consistently. 


The  imbedded  shock  is  computed  in  a  similar  way.  Let  W 
be  the  shock  velocity  and  the  values  in  front  of  the  shock  be 
denoted  by  an  index,  1;  then 
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u.  ,  =  u.-W 
Irel  1 


2  2 
fC  ,  =  (u,  ,/a.) 

rel  Irel  1 


(12) 


d  =  rM  rel-  4 


P  =  P  +  In  d/( S+1 ) 


03) 


u__,  =  4/(4+1)  u.  ,  +  a?/[(3+1)  u,_.J 


rel 


Irel 


Irel' 


3P/3W  =  -  2  u1rel/(de1) 


3U/3W  a  (M  .+1)/(S+1) 
rel 


(14) 


»  (F 1  *  (F) 

iW  -  —  [a(P  -Pl  ;)+Y(u  ,+W-u  j]/(a  3P/3W  +  y  3u/3W]  (15) 

rel 


3.  Numerical  results 


For  a  numerical  example,  a  5.56  mm  gun  was  chosen,  fc 
which  the  bullet  velocity  can  be  defined  by  the  following  equa¬ 
tions: 


b.  a  E  (t/100)5  t<t 

t  o 

bfc  s  At  +  B  +  1  /( Ct  ♦  D)  t>t 

t  o 

with  t  s  70,  E=10  and  A,B,C,D  defined  by  the  conditions: 
o 

bt  a  1.6807  ,  btt  a  .12005  at  t=70 

bfc  a  3.5  ,  btt  =  .007  at  ta125 

that  is,  Aa-0.0042,  B=4.8983,  C=-0.0145,  D=0.6739.  The  velocity 
of  the  bullet  as  a  function  of  time  is  shown  in  Fig.  1.  Some 
plots  of  pressure  distributions  are  displayed  in  Fig.  2. 

The  peculiar  shape  of  the  pressure  distribution  at  the 
highest  value  of  time  is  due  to  the  effect  of  variable  entropy 
behind  the  accelerating  imbedded  shock. 
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4.  First  phase  of  precursor  shock  evolution 


As  the  precursor  shock  reaches  the  mouth  of  the  barrel, 
we  start  the  calculation  of  the  axisymmetric  flow  around  the 
muzzle.  Typically,  the  flow  is  limited  by: 

1)  the  centerline  of  the  gun, 

2)  the  precursor  shock, 

3)  the  outer  wall  of  the  gun,  and 

4)  the  mouth  of  the  barrel. 


The  flow  is  assumed  to  be  viscous.  The  inviscid  flow 
model  would  be  both  unrealistic  and  numerically  unmanageable;  a 
preliminary  discussion  of  the  difficulties  can  be  found  in  [18]. 

A  detailed  description  of  the  equations  of  motion  is 
given  in  [17],  from  which  the  part  relevant  to  the  present  calcu¬ 
lation  is  transcribed. 


5.  Equations  of  motion 


♦ 

Let  V  be  the  velocity  vector,  t  the  time,  t  the  coeffi¬ 
cient  of  heat  conduction,  and  *  the  dissipation  term.  The  latter 
is  a  well-known  non-negative  quadratic  form  depending  on  the 
space  derivatives  of  the  velocity  components.  Its  expressions 
for  the  coordinate  system  considered  in  the  present  paper  will  be 
given  later  on.  We  will  assume  that  the  viscosity,  »,  is  a  con¬ 
stant  . 


A  Reynolds  number  and  a  Prandtl  number  can  be  defined  as 


R  -  o  _u  „x  Jv , 
e  ref  ref  ref 


P  =  c  u/k 
r  p 


(16) 


where  c  is  the  specific  heat  at  constant  pressure.  Using  the 
P 
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Equations  of  motion 


same  units  as  in  section  2,  the  equations  of  motion  for  a  viscous 
flow  (Navier-Slokes  equations)  can  be  written  in  the  form: 

Bt  *  »’-5  =  0 

DV  1  1  ii  ♦  ♦ 

=  c3  vv-v-7x7xv] 

e  (17) 


i  =  rr -ht-'). -f»20] 

e  r 


In  an  inviscid  flow,  the  terms  here  affected  by  1/R  do 

e 

not  appear.  For  a  proper  numerical  analysis  of  such  flows,  whose 
mechanics  is  governed  by  the  propagation  of  sound  waves  and  by 
the  convection  of  entropy  along  particle  paths,  it  is  convenient 
to  recast  the  equations  of  motion  into  a  form  similar  to  (4); 

DP  ♦ 

Dt  ^7-V  =° 

DV 

«■  9  VP  s  0  (18) 

=  o 

Dt 

If  the  flow  is  viscous,  the  basic  phenomena  of  wave  pro¬ 
pagation  are  still  present,  although  modified  by  the  concurrent 
effects  of  diffusion.  Therefore,  it  would  not  be  advisable  to 
drop  the  basic  integration  techniques  for  convective  terms;  we 
will  consequently  write  the  Navier-Stokes  equations  in  the  form: 

^  V  VP  *  v7  U  -  ^ 

at  ■  r  T  ,v  ’  Dt 


3  V  T  2  *  ♦ 

~  +  ^7(q  )  -  VxVxV  +  8  VP  : 
t  d 


1  ♦  - 
-5—  [rr  W.V-VxVxV] 
P  R  3 

e 


H ♦  v.vs  =  -f[(T-D*  ♦  y 

e  r 


(where  the  material  derivatives  have  been  replaced  by  partial 
derivatives,  as  costumary,  and  q  is  the  modulus  of  the  velocity). 
In  what  follows,  the  density  will  no  longer  be  used  explicitly, 
and  p  will  be  redefined  in  Eq.  (21). 
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6^  Conformal  mapping 


We  will  consider  two  types  of  flows,  both  depending  on 
two  space  variables:  two-dimensional  (plane)  flows  and  axi- 
symmetric  flows.  Two  Cartesian  coordinates,  x  and  y,  will  be 
used  for  two-dimensional  flows.  The  same  symbols  will  be  used 
for  the  axial  and  radial  coordinate,  respectively,  in  any  meri¬ 
dional  plane  of  an  axi-symmetric  flow.  The  (x,y)-plane  will  be 
called  the  physical  plane.  We  must  provide  a  computational  grid 
with  a  family  of  lines  running  from  the  centerline  to  the  wall 
and  another  family  of  lines  running  from  the  mouth  to  the  3hock. 
To  this  effect,  we  begin  by  introducing  a  complex  variable 

z  =  x  +  i  y  (20) 

and  assume  that,  in  general,  the  physical  plane  will  be  confor¬ 
mally  mapped  onto  an  auxiliary  plane,  described  in  terms  of  a 
complex  variable  C, 


c 


ie 


P  e 


(21) 


Such  a  mapping  is: 


z  =  (rQ/ir)[(z^-  1/z^)/2  -  log  z^  -  i*] 
z  +  1/z1  =  2B  (C  +  1/C) 


(22) 


where  rQ  is  the  outer  radius  of  the  barrel,  and  B  is  determined 
to  assure  correspondence  between  C=1  and  z  =  -i.  Lines  of  con¬ 
stant  p  and  of  constant  9  in  the  C-plane  are  mapped  onto  the  grid 
3hown  in  Fig.  3.  If  a  9  =  constant  line,  different  from  BD  (for 
example,  FG),  is  used  as  a  boundary,  the  flow  in  the  physical 
plane  will  occur  in  the  divergent  duct  having  FG  as  a  wall.  Let 


g 


•  S  ■  o*1* 


(23) 
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If  I  and  0  are  the  unit  vectors  parallel  to  the  x-  and  y-axis, 
respectively,  and 


V  =  U  1  ♦  V  3 


(32) 


it  follows  that 


1.1=2  ,  £.0=2  ,  3.I=-S  ,  3.0=2  (33) 

U=u£-v3,  V=uS+v2,  u=U2+V2,  v=-U2+V2  (3*0 

For  any  element  ds  jji  space,  we  have 

.  2  1,2  £>2  Jfl2  .  2,  2 

ds  =  —gd p  +  d9  ♦  jy  dx^ 

G  G  (35) 

♦ 


if  x^  is  the  angular  coordinate  in  an  axi-symmetric  problem  and  j 
is  a  multiplier,  equal  to  0  for  two-dimensional  flows  and  to  1 
for  axi-symmetric  flows.  Consequently,  the  basic  vector  opera¬ 
tors  in  (19)  can  be  expressed  as  follows: 

VP  =  G(Fp£  +  j  Pej)  (36) 

T"?.  -  (§V  (37) 

2 

VxVxV  =  ^-[(^-^  -  (§)0](v£-u3)  (38) 


v.v  =  5-[(£^) 

P  G  p 


<5>e> 


(39) 


and 


*ft<eire22>2*<'22-'33>2*<e33-'n>2)  (“0) 


with 


*11  *  G(UP  +  p  V 


e 22  =  flV  u(Ml)]  •  e33  =  j  y 
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el2  =  J(u9-^2)] 


Computational  plane 


7.  Computational  plane 


The  C-plane,  however,  is  not  the  computational  plane. 
Let  o=c(9,t)  be  the  image  of  the  precursor  shock  in  the  C-plane. 
A  normalized  computational  grid  is  then  obtained  by  letting: 


X  =  (0  -  1)/(c  -  1) 


8  s 


1  9  ,tanh[a(Y-1 )  ] 

2  o'  tanh  a  +  2 


T  =  t 


(48) 


where  9q  is  the  constant  value  of  9  which  defines  the  rigid  wall 
and  which,  in  the  present  case,  should  be  equal  to  zero.  The 
computation  is  performed  in  the  (X,  Y)  plane,  using  a  rectangular 
grid  with  evenly  spaced  lines;  it  must  be  noted,  however,  that 
X=constant  lines  and  Ysconstant  lines  on  the  physical  plane  are 
generally  not  orthogonal  to  each  other,  since  the  boundary 
0=c(8,t)  is  not  constant  with  respect  to  8.  The  lack  of  ortho¬ 
gonality,  however,  does  not  affect  the  accuracy  of  the  results. 

Since  here  two  boundaries  are  defined  by  constant  values 
of  8  and  of  the  other  two  boundaries  one  is  defined  by  a  constant 
value  of  0  but  the  other  by  a  value  of  o  which  depends  on  8  and 
t,  it  will  be  necessary  to  define  X  not  only  as  a  function  of  0 
but  of  8  and  t;  therefore,  we  will  write: 

X  =  X( P , 9 , t ) 

Y  =  Y ( 9 )  (49) 


T  =  t 


Letting 


Computational  plane 


and 


b12  =  yGY0/p  ,  g12  =  yGX9/o  ,  a21  =  Gexp 


h3  =  G9X0/P  ,  b21  =  G6Y0/P 


F  =  C (Y-l )♦  +  )-V26]/(pRe) 


c,  =  yE  -  F 

If  ‘1VV°  -  j  f ' 11 


c3  ■  -uD  -  |f  [fe/c*a.  *  J  It01 


the  equations  of  motion  are  recast  in  the  form: 


(50) 


(51) 


(52) 


PT+a11PX+b11-Y 


P,.+a  ,0uv+b,  ,,vv-*-g 


12  X  12  Y 


i2vx+cr° 


UT+a11UX+b11UY+a21PX+C2=° 


Va1lVb1lVb2lVh3PX*V° 


ST+al1SX+t>1  1SY  =  F 


(53) 


It  is  easy  to  see  that  the  first  three  equations  (53)  and  the 
first  three  equations  (65)  have  the  same  form  as  Equations  (23) 
of  [3].  The  integration  procedure  explained  in  [3]  can  be  ap¬ 
plied.  We  expect  the  integration  technique  to  provide  a  very 
good  estimate  of  convection  and  wave-propagation  effects  which 
are  common  to  inviscid  and  viscous  flows. 
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8.  Rigid  boundaries 


For  a  viscous  flow,  the  velocity  at  a  rigid  boundary  is 
assumed  to  vanish.  There  are  no  difficulties,  thus,  in  the 
determination  of  u  and  v.  Pressure  and  entropy  require  more 
care.  We  consider  here  the  case  of  an  isothermal  wall;  the  tem¬ 
perature,  9  is  a  prescribed  constant.  Because  of  (17),  P  and  S 
are  linearly  related;  therefore, 


DS 

DT 


(  1  — Y  ) 


DP 

DT 


and  consequently,  the  first  of  (46)  must  be  replaced  by 

P„  ♦  G(uP  ♦  P .)  +  G(u  +  ;  v.)  =  0  (54) 

t  0  p  9  POO 


having  taken  into  account  that  E  vanishes  identically.  The  equa¬ 
tions  to  be  solved  are  still  (53),  provided  that 

a12’  b12’  g1l  and  812  are  div*ded  by  Y  and  ci  is  set  ecIual  to 
zero . 


Once  the  pressure  has  been  determined,  S  is  obtained  from 
(17).  Therefore,  S  is  not  calculated  directly  as  a  result  of 
dissipation  and  heat  transfer  in  the  flow,  because  the  condition 
of  constant  wall  temperature  implies  some  external  action,  and 
(17)  gives  us  the  final  outcome  of  such  action  and  the  local 
variation  of  pressure. 


9.  Integration  procedure 


The  equations  of  motion  (53)  are  integrated  following  the 
general  guidelines  of  Section  6  in  [33.  We  define 

C1=811UY  +  C1'  C2=bnUY  +  h2PY  +  C2 
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C1=g12VX 


(55) 


c2'a1 1  vx  ■"  h3Px 


'3 


X  Y 

After  finding  the  characteristic  slopes,  X^,  X ^  (i=1,2)  from 


.X 

a11  i 

a21 

O 

1 

H-  *-< 

b21 

a „  „ 

a  -xX 

=0  . 

h  h  -XY 

12 

11  i 

12 

11  i 

=0 


(56) 


and  letting 


.X 

A  _  i 

1  '  xX_xX 

2  1 


b  -XY 

Bi  =  -r~r 

X2~X\ 


(57) 


DX. 

1J 


ILL  dy 
x  ,x  •  Uij 


_LL 


X  -X 

2  1 


Y  Y 
X  -X 
2  1 


the  equations  to  be  integrated  are: 


PT  *  VfofVs'M 


UT  *  D21(X2PX2-11PX1)  *  A1X2UX2-A2X 5UX1 


t>1!2<X2“x2-il!uX.>  '  C1  =° 

+  cx  =0 


(58) 


and 


P?  *  B.'IPY,-B2'2PY2  *  “V^YJ^YI1  *  *° 

VT  *  D21<l2Pr2-*'[Prt>  *  B1X2VY2-B2X 1VY1  *  C2  =° 


(59) 


with 


X  Y 

p  _  pA  +  p; 

T  r  rT 


(60) 


For  discretization,  the  space  derivatives  are  classified 
into  three  categories:  (i)  the  one3  explicitly  appearing  in  (58) 
and  (59),  (ii)  the  ones  explicitly  appearing  in  (55),  and  (iii) 
the  ones  appearing  in  (51)  and  (52).  The  derivatives  of  the 
first  category  are  discretized  according  to  the  rules  (1*0  and 
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(15)  of  [31;  the  derivatives  of  the  second  category  according  to 
the  rules  (31*)  and  (35)  of  [33;  an  the  derivatives  of  the  third 
category  are  approximated  by  ordinary  centered  differences,  be¬ 
cause  the  physical  nature  of  the  terms  which  they  affect  is  dif¬ 
fusion.  Few  exceptions  to  the  general  rules  are  necessary  for 
points  on  rigid  boundaries  or  next  to  rigid  boundaries.  For 
points  on  rigid  boundaries,  the  alternate  two-point-three-point 
approximation  is  always  taken  using  points  inside  the  flow  field. 
For  points  next  to  rigid  boundaries,  if  use  of  a  point  located 
behind  the  wall  were  required  in  a  three-point  approximation,  the 
latter  is  substituted  by  a  two-point  formula. 


10.  Numerical  treatment  of  shocks 


Shocks  are  generally  present  in  a  compressible  flow 
field,  either  as  boundaries  or  imbedded  in  the  flow.  In  both 
cases,  they  can  be  treated  numerically  in  the  general  framework 
of  the  computational  technique  described  in  the  previous  Sec¬ 
tions.  We  begin  with  some  general  considerations. 

Let  f)  and  t  be  the  unit  vector  normal  and  tangential  to  a 
shock  at  any  of  its  points,  Q,  respectively: 

R  =  Nt£  ♦  N2j  ,  ?  =  -N2£  +  N1 j  (61) 

♦ 

W  the  shock  velocity. 


W  =  W  R 


(62) 


and  u,  v  the  velocity  components  along  ft  and  ?,  respectively: 
u  =  uN1+vN2  u  s  GNrvN2 

v  =  -u  N2+vN  ^  v  =  uN^+vN^ 

The  N-component  of  the  velocity  relative  to  the  shock  is 


(63) 
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urel=u-W  (64) 

The  relative  normal  Mach  number  on  the  low  pressure  side  of  the 
shock  is 


nlrel  "  Y0, 


The  Rankine-Hugoniot  conditions  are: 


?2  =  P,  +  In 


2yM  ,-y+1 
n  Irel 


Y@ 

Y— 1  -  2  T  1 

u2rel  “  Y+1  U1rel  +  Y  +  1  u,  , 

Irel 


We  will  need  the  derivatives  of  P 2  and  u2j-ei  with  respect 


Su^i-df-De, 


Izi  ^  T91 

Y+1  +  Y+1  ~2 

U  Irel 


Let  us  assume  that  the  shock  is  oriented  in  the  general 
direction  of  the  Xsconstant  lines;  therefore,  it  can  be  defined 
by  its  intersections,  Xg ,  with  Y=constant  lines.  At  each  inter¬ 
section,  we  consider  two  points,  one  on  the  low-pressure  3ide  and 
the  other  on  the  high-pressure  side.  The  point  on  the  low- 
pressure  side  must  be  updated  by  using  information  proceeding 
from  that  side  only.  The  point  on  the  high-pressure  side  must  be 
updated  by  using  information  from  both  sides.  The  information 
proceeding  from  the  low-pressure  side  must  satisfy  the  Rankine- 
Hugoniot  conditions;  the  information  proceeding  from  the  high- 
pressure  side  is  carried  to  the  shock  along  a  characteristic. 
The  acceleration  of  the  shock  results  from  the  compatibility  of 
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the  different  information.  To  obtain  the  characteristic  equa¬ 
tion,  we  should  rewrite  the  equations  of  motion  in  a  frame  rela¬ 
tive  to  the  moving  shock,  where  the  Y-coordinate  of  the  shock  is 
unchanged  but  its  X-coordinate  follows  the  shock  in  its  motion. 
Therefore,  we  will  introduce  a  new  set  of  coordinates. 


X  =  X  -  X  (Y.T) 


e  =  Y 


t  =  T 


and  rewrite  (65)  in  the  form: 


PT+a11PX+b11Pc+ai2UX+b12Vc+Y12Vx+C1=0 


UT+ai1UX+bnUc+a21PX+h2Pc+C2=° 


Vx+ai1Vx"b11Ve+b21Pc"n3PX+C3=° 


S  +a, ,S  +b, ,S  =F 
t  11  x  He 


where 


®irairb1 1XsY'XsT  ’  al2=a12-g1 1XsY  •  a21'a2rh2XsY 


Y12=g12"b12XsY  *  n3=h3_b21X3Y 


A  characteristic  equation  may  now  be  obtained,  using  x  and  t  as 
basic  independent  variables: 

(<*11-A)(PT-*-XPx)-a12(uT+Aux)-x12(vT+Avj()*R  =0  (7D 

where  X  is  the  solution  of  the  equation 


°irx  a2i  n3 


°i2  airA  0  V  0 


0  «n-A| 
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that  is. 


X  =  air(al2a21+Y12n3) 

and  R  comprises  all  the  remaining  terms. 


1/2 


(73) 


Note  now  that 


“12  =  YGvN1  ,  y12  =  yGvN2 


421  ~  G9vNi  •  n3  -  G8vN2 


where 


Therefore , 


yG 


a,2u  +  Y,2V  =  YGvu 


“12ut  ♦  1f12VT  s  1rGvuT-'rGv(uNlT-t-vN2T) 


“12UX  +  T12VX  YGVUX 


Instead  of  (71)  we  can  write: 


where 


(«1  j-X  )(PT+*P  )-yGv  (ut+Xux  )+R  1  =  0 


R i  =  R+yGv (uN1t~vN2t ) 


A  further  simplification  is  obtained  by  observing  that 

1  /2 

<*,  -X  (a  a  +y  n,) 

1 1  _  +  12  21  12  3  _  +a 

yGv  "  “  yGv  ~y 


(74) 


(75) 


(76) 


(77) 


(78) 


because  of  (74)  and  (75),  so  that  (77)  can  be  written  in  the 
form : 
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(PT+xPx)^T>xnx+Ri  =  0  (79) 

The  values  of  ?T  and  ut  ,  obtained  by  integrating  the 
Navier-Stokes  equations  (69)  using  information  from  the  high- 
pressure  side  of  the  shock  only,  must  satisfy  (79).  On  the  other 
hand,  the  exact  solution  of  the  flow  problem  in  the  presence  of 
the  shock,  which  accounts  also  for  the  information  from  the  low- 
pressure  side  and  the  Rankine-Hugoniot  conditions,  must  satisfy 
(79)  as  well.  Therefore,  if  we  write  (79)  twice,  once  for  the 
t -derivatives  as  obtained  from  the  Navier-Stokes  equations  (that 
is,  using  for  the  shock  point  the  same  integration  procedure 
which  is  applied  to  other  points)  and  again  for  the  exact 
t -derivatives  and  we  subtract  one  equation  from  the  other,  the 
simple  result  is  obtained: 


~  -(NS) 

ut-uT 


=  0 


(80) 


where  the  derivatives  obtained  from  the  Navier-Stokes  equations 
are  labeled  (NS)  and  the  exact  derivatives  are  unlabeled. 


The  latter  derivatives  can  obviously  be  expressed  in  the 

form : 


r 


i 


3P 

P  =  P»  -rrr—  W 

T  T  3W  T 


3u2 

UT  =  “?  +  WT 


(81) 


where  f*  is  a  derivative  computed  considering  W  as  a  constant. 
The  acceleration  of  the  shock  is  thus  obtained: 


+  ,D(NS)  .-(NS) 

-a(PT  -P*)+y(ut  -u»] 

±a  3P2/3W  +  y  3u2/3W 


(82) 


Since  both  the  starred  derivatives  and  the  derivatives  indicated 
by  (NS)  are  computed  using  the  same  initial  values,  (82)  can  be 
replaced  by 


wrat  = 


.  ,D(NS)  D#,  .-(NS) 

ta(P  — P  * ) +y ( u  -u*! 

-a  3P2/3W  +  y  3u2/3W 


(83) 
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where  P*  and  u*  are  the  values  on  the  high-pressure  side  of  the 
shock  obtained  from  updated  values  on  the  low-pressure  side  by 
applying  the  Fankine-Hugoniot  conditions  to  a  shock  whose 
geometry  has  been  updated  but  whose  velocities,  W,  are  still  the 
same  as  at  the  beginning  of  the  integration  step.  This  result 
[4]  is  remarkable  since  it  provides  an  extremely  simple  method 
for  calculating  the  shock  acceleration  although  it  relies  on  the 
same  basic  concepts  which  have  been  shown  to  be  necessary  for  a 
proper,  physically  correct,  handling  of  shocks  [5]. 

We  will  now  consider  first  the  case  where  the  shock  is  a 
boundary,  moving  into  a  gas  at  rest,  and  then  the  case  of  an  im¬ 
bedded  shock. 

In  the  case  of  a  shock  moving  into  a  gas  at  rest,  let  us 
assume  that  the  shock  is  a  right  boundary  of  the  computational 
field.  It  is,  thu3,  defined  by  X=1;  in  this  case,  x,  e  and  t 
coincide  with  X,  Y  and  T,  respectively  and  a^ra^,  a21=a2r 

ai2=a12'  Y12=®12'  n3=^3‘  ckaracteristic  reaching  the  shock 

from  the  high-pressure  side  is  a  right-running  characteristic , 
and  in  the  preceding  equations,  whereas  a  ±  appears,  the  upper 

sign  must  be  used.  Note  in  addition  that  u,=0  and  v,=0;  there- 

2  2  11 
fore,  u.  . =-W  and  M  ,  rW  /y .  In  this  case,  obviously,  the 
lrel  n  lrei 

low-pressure  side  values  are  known  without  any  need  for  comput¬ 
ing;  the  (NS)-values  on  the  high-pressure  side  are  obtained  to¬ 
gether  with  and  using  the  same  procedure  as  for  interior  grid 
points . 


For  an  imbedded  shock,  whose  location  does  not  generally 
coincide  with  a  grid  line,  we  use  a  simplified  procedure  to  ob¬ 
tain  the  values  at  the  shock  on  the  low-pressure  side  and  the 
(NS)  values  on  the  high-pressure  side.  On  the  low-pressure  side, 
instead  of  integrating  (69)  (which  would  require  a  special,  and 
not  easy,  redefinition  of  approximations  for  the  x-  and 
c-der ivatives) ,  we  simply  assume  that  the  values  at  the  shock  can 
be  extrapolated  from  the  two  adjacent  grid  points  on  the  same 
y=constant  line,  both  at  the  end  of  the  predictor  and  the  correc¬ 
tor  level.  On  the  high-pressure  side,  we  assume  that  the  T- 
derivative3  on  the  shock  are  equal  to  the  T-der ivatives  at  the 
grid  point  next  to  the  3hock  on  the  same  Yrconstant  line.  Note, 
however,  that 
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f  =  f 
t  1 T 


f  X  _ 
X  sT 


(84) 


for  any  function  f.  The  values  of  f^.  at  the  shock  are  assumed  to 
be  the  same  as  at  the  next  grid  point  on  the  high-pressure  side, 
on  the  same  Y=constant  line.  The  values  of  f  are  approximated 
as  follows: 


f  =[f.-f 
X  A  s 


(85) 


where  e=(X^-Xg)/4X  and  A  is  the  grid  point  next  to  the  shock  on 
the  high-pressure  side  and  B  is  the  next  grid  point.  This  formu¬ 
la  provides  a  smooth  transition  when  the  shock  crosses  an 
X=constant  line. 


For  an  imbedded  shock,  thus,  the  calculation  proceeds  as 
follows:  In  the  predictor  stage,  after  updating  all  grid  points, 
the  low-pressure  side  of  the  shock  is  obtained  by  extrapolation 
and  (84)  is  also  applied  to  P,  u,  v  and  S.  The  values  on  the 


high-pressure  side  are  updated  by  adding 


f/t 


to 


the  initial 
( NS ) 

values  of  f;  the  values  so  obtained  are  the  predicted  i  .  The 
predicted  f*  are  obtained  by  applying  (66)  to  the  predicted 
values  on  the  low-pressure  side.  Then,  (83)  i3  applied  and  W  is 
temporarily  updated,  but  its  original  value  is  retained  in 
storage.  The  geometry  of  the  shock  is  updated,  considering  that, 
in  virtue  of  (62),  (28)  and  (33), 


Pgt  =  G  W/N1  (86) 

and  using  the  approximations: 

»,<'*«)  =  »,<*>  *  p«“  *  <87) 

where  the  second  derivatives  are  obtained  by  differentiating 
(86).  In  the  corrector  stage,  the  procedure  outlined  above  for 
the  predictor  stage  is  repeated  through  the  application  of  (84). 
The  updating  of  the  values  on  the  high-pressure  side  is  obtained 
by  adding 


2 


(f. 


f ( pred ) j 


to  the  predicted  values.  The  corrected  f*  are  obtained  by  apply- 
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ing  (66)  to  the  corrected  values  on  the  low-pressure  side,  with 
the  original  value  of  W.  Then,  (83)  is  applied  again  and  W  is 
definitively  updated.  The  Rankine-Hugoniot  conditions  (66)  are 
applied  once  more  to  the  corrected  values  on  the  low-pressure 
side  using  the  final  value  of  W,  to  obtain  the  final  values  on 
the  high-pressure  side.  The  entropy  is  also  computed  from 


S-  =  S.  +  P»-P.-  Yln(u.  ,/u_  ,) 

2  1  21  Irel  2rel 


(88) 


At  this  stage,  it  is  convenient  to  correct  the  values  at 
the  grid  point  next  to  the  shock  on  the  high-pressure  side. 
Pressure  and  velocity  components  are  interpolated  from  the  values 
at  the  shock  and  the  following  grid  point.  Entropy  is  also  in¬ 
terpolated  considering  that  it  is  carried  along  a  streamline. 


11.  Details  of  the  calculation 


The  calculation  is  started  at  a  small,  positive  value  of 
t,  when  the  precursor  shock  has  already  moved  from  the  muzzle 
mouth.  Initially,  the  shock  is  assumed  to  lie  on  a  P=constant 
line  (with  the  value  of  the  constant,  cq,  slightly  greater  than 
1),  and  all  parameters  pertinent  to  the  shock  are  assumed  equal 
to  their  values  at  t=0.  Note  that,  since  the  flow  behind  the 
3hock  is  uniform,  the  shock  Mach  number  defines  pressure,  veloci¬ 
ty  and  entropy.  Uniform  flow  is  assumed  between  the  muzzle  mouth 
and  the  shock,  if  the  gas  is  inviscid.  The  velocity  component,  v 
is  made  equal  to  zero  throughout. 

Initially,  the  computational  region  is  divided  only  into 
two  strips  along  the  X-axis.  Along  the  Y-axis  we  consider  as 
many  partitions  as  necessary  to  provide  sufficient  resolution. 

Since  the  flow  is  viscous,  we  need  some  modifications  to 
the  initial  conditions  near  the  wall,  to  account  for  the  vanish¬ 
ing  of  the  velocity  at  the  wall.  The  precursor  shock  cannot  reach 
the  wall;  the  perturbation  front,  elsewhere  in  the  form  of  a 
shock,  becomes  a  characteristic  at  the  wall.  Therefore,  on  the 
initial  o=c  line  which  represents  the  shock,  P  is  assumed  equal 
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to  zero  at  the  wall.  On  the  next  wall  point,  P  is  taken  equal  to 
one  half  of  its  value  behind  the  shock.  The  wall  temperature  Is 
assumed  equal  to  1  at  any  time.  The  entropy  is  defined  accord¬ 
ingly.  The  u-velocity  component  is  taken  equal  to  one  half  of 
its  value  behind  the  shock  along  the  second  9=constant  line  from 
the  wall. 

From  the  initial  conditions,  the  calculation  proceeds  as 
detailed  above.  We  note  that  the  centerline  is  computed  as  a 
symmetry  line,  not  as  a  rigid  wall.  For  computational  purposes, 
the  shock  geometry  is  prolonged  to  the  wall  with  a  constant  value 
of  c,  and  the  wall  point  is  computed  as  any  other  wall  point, 
taking  advantage  of  the  fact  that  the  state  of  the  gas  in  front 
of  the  perturbed  region  is  known.  A  similar  procedure  is  au¬ 
tomatically  applied  to  any  other  point  on  the  perturbation  front, 
if  the  3hock  happens  to  lose  its  strength  completely:  that  is, 
the  point  on  the  perturbation  front  is  computed  as  any  interior, 
shockless  point,  taking  advantage  of  the  fact  that  the  state  of 
the  gas  in  front  of  the  perturbed  region  is  known. 

As  the  calculation  proceeds,  the  number  of  grid  intervals 
in  the  X-direction  is  doubled  every  time  ( c— 1 )  on  the  centerline 
exceeds  1.4  times  its  initial  value  or  the  value  it  had  at  the 
previous  doubling,  until  the  total  number  of  intervals  is  16. 

Certain  features  of  the  flow  are  common  to  the  inviscid 
and  viscous  models.  An  expansion  appears  from  the  beginning  near 
the  inlet  of  the  duct;  the  region  of  maximum  expansion  moves  from 
left  to  right,  but  at  a  slower  speed  than  the  precursor  shock. 
Consequently,  even  with  the  precursor  shock  losing  strength,  the 
pressure  behind  it  remains  higher  than  the  lowest  pressure  al¬ 
ready  attained  along  the  duct.  The  particles  are  accelerated  and 
then  decelerated  again.  A  compression  wave  appears,  which  tends 
to  steepen  up,  as  every  compression  wave  does,  and  another  shock 
results  eventually. 

In  the  presence  of  viscosity,  the  recompression  wave  and 
the  secondary  shock  strongly  interact  with  the  boundary  layer  in 
the  process  of  formation.  The  latter  thickens  and  separates  very 
soon.  Between  the  main  stream  and  the  wall,  a  wide  dead-water 
region  appears,  where  the  pressure  tends  to  equalize  the  ambient 
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pressure  (in  front  of  the  precursor  shock).  The  separated  flow 
never  reattaches  and  becomes  a  plume.  From  the  separation  point 
on,  the  plume  is  insensitive  to  the  wall  geometry  and  the  flow 
inside  it  is  essentially  inviscid.  The  results  of  the  calcula¬ 
tion  of  a  steady,  inviscid  flow  in  a  plume  can  be  used  to  judge 
whether  the  viscous  calculation  approaches  its  theoretical  asymp¬ 
tote,  and  how  well. 

Three  cases  have  been  calculated,  with  different  values 

of  9^  and  of  the  initial  conditions.  The  first  deals  with  a  high 

value  of  9  by  virtue  of  which  the  flow  occurs  in  a  divergent 

(two-dimensional)  channel.  The  second  and  the  third  use  9  =0. 

o 

The  second,  however,  does  not  contemplate  initial  conditions  con¬ 
sistent  with  the  flow  in  the  barrel  described  in  Section  3;  rath¬ 
er,  it  has  initial  conditions  consistent  with  a  shock-tube  exper¬ 
iment  [19].  The  third  case  is  consistent  with  the  barrel  flow 
described  in  Section  3. 


12.  First  example 


In  the  first  example,  9q  has  been  taken  equal  to  1.2. 

and  rQ  equal  to  2.2603.  Therefore,  the  calculation  attempts  to 

describe  the  flow  in  a  diverging  channel,  inserted  at  the  end  of 

a  straight  shock  tube  (Fig.  4).  The  flow  is  assumed  to  be  two- 

dimensional.  The  stretching  parameter,  a,  has  been  taken  equal 

to  2,  and  the  sector  between  9=9  and  9=n/2  is  divided  into  16 

o 

intervals.  Consequently,  we  obtain  a  fair  accumulation  of  grid 
lines  near  the  wall,  and  still  work  with  a  reasonably  small 
number  of  lines.  It  is  clear,  however,  that  the  resolution  is 
well  below  the  limits  which  are  usually  recommended  for  a  good 
description  of  Reynolds  number  effects;  the  lack  of  an  adequate 
resolution  is  particularly  felt  in  the  vicinity  of  the  plume 
shear  layer.  Nevertheless,  the  present  results  are  very  en¬ 
couraging,  just  because  very  good  qualitative  results  are  ob¬ 
tained  with  such  a  coarse  mesh.  The  Reynolds  number  used  in  this 
case  is  extremely  high  (10  ),  so  that  the  flow  should  be  con¬ 
sidered  as  inviscid  practically  everywhere,  except  across  the 
boundary  layer.  We  know,  however,  as  we  3aid  in  the  preceding 
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Fig.  ** 


section,  that  the  steady  configuration  of  the  flow  in  a  channel 
whose  area  increases  beyond  any  limit  is  a  plume,  that  is,  a 
separated  flow.  The  plume  exists,  regardless  of  the  value  of  the 
Reynolds  number;  only  the  thickness  of  the  shear  layer  limiting 
the  plume  reveals  what  the  Reynolds  number  is.  In  the  present 
case,  we  should  expect  the  formation  of  a  plume  with  a  limiting 
shear  layer  of  practically  vanishingly  small  thickness;  but  we 
also  know  that  such  a  picture  is  going  to  be  distorted  by  the  ar¬ 
tificial  Reynolds  number  of  the  computational  me3h,  which  is 
smaller  by  orders  of  magnitude. 

The  Prandtl  number  is  taken  equal  to  1,  and  r  equal  to 
l.1*.  The  temperature  at  the  wall  is  assumed  to  be  equal  to  1, 
that  is,  to  the  temperature  of  the  gas  at  rest.  The  value  of  ? 
at  the  inlet  i3  i.6«87.  Consistent  values  of  u  and  S  are 
usl.6551,  SsO. 1702.  The  Mach  number  at  the  inlet  is  barely 
supersonic,  Msl.juo. 

The  progression  of  the  precursor  3hock  and  the  flow  evo¬ 
lution  behind  it  are  shown  in  Fig.  5  by  sets  of  isobars  drawn  at 
different  instants  of  time.  The  figures  represent  the  lower  half 
of  the  channel.  Cn  the  symmetry  line,  notches  denote  intervals 
of  unit  length.  The  isobars  correspond  to  constant  values  of  P, 
spaced  as  indicated  oy  the  value  ;f  OREF  m  each  picture;  the 
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lines  marked  with  0  are  sonic  lines. 

At  t=0.6665,  the  expansion  has  barely  begun  and  is  still 
very  mild;  at  t=1.3718  it  is  stronger  (the  lowest  value  of  P  is 
less  than  0.5)  and  the  recompression  is  clearly  taking  shape. 
Within  the  first  unit  length  of  the  channel,  the  interaction  of 
up-  and  down-running  characteristics  is  very  clear  (an  obvious 
consequence  of  M  being  very  close  to  1  at  the  start);  the 
recompres3ion ,  therefore,  appears  almost  as  unidimensional.  At 
t=2.  1100,  the  recompression  is  clearly  piling  up  at  two  unit 
lengths  from  the  entrance;  the  lowest  value  of  P  i3  below  the 
ambient  value.  As  indicated  by  the  modified  shape  of  the  3onic 
line,  a  small  separated  region  appears;  this  effect  is  much 
clearer  in  Fig.  6,  where  pictures  of  streamlines  are  shown.  A 
shock  is  soon  formed,  and  it  is  evident  in  the  isobar  plot  at 
t  =  3.i*77I*  (the  thickness  of  different  arcs  of  the  shock  is  propor¬ 
tional  to  its  local  strength).  The  second  drawing  in  Fig.  6 
clearly  shows  the  separation  of  the  flow  at  t=4.091^.  The 
streamlines  in  the  separated  region  should  not  be  interpreted  as 
picturing  a  flow  as  strong  as  the  main  flow;  as  a  matter  of  fact, 
the  velocities  in  the  separated  region  are  extremely  small,  in 
general . 


The  isobar  picture  at  t=4.5097  3hows  a  changing  pattern; 
the  recompression  3hock  shrinks  towards  the  symmetry  line,  and 
all  i3obar3  with  negative  values  of  P  (that  is,  below  ambient 
pressure) ,  as  well  as  the  P=0  line,  show  a  definite  tendency  to 
bend  horizontally  in  the  direction  of  the  flow.  At  t=6.2530,  the 
qualitative  description  of  the  plane  is  correct,  but  resolution 
ha3  become  so  poor  that  most  of  the  quantitative  results  are  dis¬ 
torted  . 


Nevertheless,  it  may  be  interesting  to  see  how  close  the 
flow  is  to  a  steady,  invi3cid  plane,  at  this  early  stage  of  evo¬ 
lution.  To  this  effect,  a  program  was  written  to  describe  a 
steady,  two-dimensional,  shockless  inviscid  plume  evolving  from 
the  same  inlet  conditions  along  the  same  channel  geometry  until 
the  ambient  pressure  is  reached,  and  then  continuing  inside  the 
region  delimited  by  a  psconstant  streamline.  The  corresponding 
isobar  pattern,  with  values  of  P  spaced  0.25  apart,  is  shown  in 
the  upper  part  of  Fig.  7.  Isobars  from  the  last  plot  of  Fig.  5, 
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Fig.  6(b) 


interpolated  at  the  same  values  of  P,  are  shown  in  the  lower  part 
of  Fig.  7.  At  this  stage,  with  no  supporting  calculations  on  a 
finer  grid,  it  is  not  possible  to  state  whether  the  overexpansion 
shown  by  the  latter  is  a  temporary  effect,  tending  to  disappear 
in  a  further  stage  of  evolution,  or  simply  the  result  of  lack  of 
resolution  (note  that  the  steady  plume  calculation  ha3  2C  inter¬ 
vals  along  any  line  orthogonal  to  the  symmetry  line  and  hundreds 
of  steps  along  the  channel,  whereas  the  lower  figure  has  been 
computed  on  the  grid  of  Fig.  ^ ) . 


22*  Second  example 

A  second  example  was  run,  based  on  experiments  performed 
on  a  muzzle  jet  flow  simulator  [’9^.  The  experimental  device 
consists  of  a  shocK  tube,  opening  into  a  low-pressure  chamber. 
The  flow  i3  ax isymmetr ical .  The  outer  radius  of  One  shock  tube 
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la  2.08  times  the  inner  radius.  The  latter  is  19  mm.  At  sea 
level,  the  speed  of  sound  is  360  m/sec;  therefore  our  reference 
velocity  is  304.25  m/sec.  Assuming  that  the  kinematic  viscosity 
of  air  is  1.6x1Q-^  m^/sec,  the  Reynolds  number  is  approximately 
36CO.  We  take  the  Prandtl  number  equal  to  1  and  y  =1.4.  We  will 
try  to  simulate  the  experiment  shown  in  Fig.  7  of  [19];  there¬ 
fore,  we  will  use  a  value  of  P  at  the  shock  tube  mouth  equal  to 
3.40.  Consistent  values  of  u  and  S  are  u=4.8250,  3  =  1. '387.  The 
tfach  number  is  1.6680. 

The  evolution  of  the  flow  is  shown  in  Fig.  3  by  sets  of 
isobars  drawn  at  different  instants  of  time,  as  it  was  made  in 
Fig.  5.  Our  results  can  be  compared  to  Fig.  7  of  [19],  taking 
into  account  that  our  reference  time  is  0.0624  msec. 

In  Fig.  9,  the  location  of  the  precursor  shock  and  of  the 
recompression  shock  on  the  centerline  is  plotted  as  a  function  of 
time;  experimental  results  from  Fig.  14  of  [ '?]  are  also  plotted 
in  the  same  figure. 
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14 .  Third  example 


Finally,  a  calculation  was  made  to  simulate  experiments 
by  Schmidt  and  Shear  on  the  M-16  rifle  [20].  Assuming  the  inner 
radius  of  the  gun  as  the  unit  length,  the  outer  radius  is  2.2603 
and  the  distance  to  be  traversed  by  the  no3e  of  the  bullet  before 
it  reaches  the  muzzle  is  150.82.  We  are  now  using  a  program 
which  couples  the  inner  and  outer  flow  calculations,  and  we 
determine  the  inner  flow  as  in  the  first  example.  The  conditions 
outside  the  gun  at  the  instant  of  firing  are  assumed  as  the  stan¬ 
dard  sea  level  conditions.  The  pressure  ratio  immediately  after 
the  precursor  shock  at  the  instant  of  its  exit  is  11.61  (slightly 
less  than  the  one  mentioned  in  [20];  the  difference  is  due  to  the 
fact  that  the  flow  velocity  behind  the  shock  when  it  reaches  the 
mouth  of  the  barrel  is  lower  than  the  bullet  velocity  when  the 
base  of  the  bullet  separates  from  the  muzzle).  To  compare  calcu¬ 
lations  and  experiments,  we  note  first  that  our  reference  time 
[equal  to  the  reference  length,  2.7813  mm.  multiplied  by  the 
square  root  of  y  and  divided  by  the  speed  of  sound  of  the  gas  at 
rest  (360  m/sec)],  is  9.1  usee.  The  precursor  shock  exits  from 
the  barrel  at  t=108,  that  is,  .983  msec  after  firing.  To  reach 
that  stage,  about  700  computational  steps  were  necessary,  in 
which  the  one-dimensional  internal  flow  alone  was  computed . 

A  set  of  isobar  plots,  at  different  instants  of  time,  is 
presented  in  Fig.  10.  The  structure  of  the  precursor  blast  com¬ 
pares  well  with  the  experimental  evidence.  The  plot  at 
t  = 1 1 0 .  1 26 3  and  the  one  at  t=113.3647  are  very  close  to  the  pic¬ 
tures  3hown  as  a)  and  b)  in  Fig.  3  of  [20],  which  correspond  to 
the  3ame  instants  of  time.  A  striking  agreement  is  shown  by  the 
locations  of  blast  shock  and  recorapression  shock  on  the  center- 
line  (Fig.  12).  The  straight  lines  in  the  figure  are  the  ones 
defined  in  [20]  as  the  best  fit  for  their  experimental  values. 
Note  that  the  instant  at  which  the  base  of  the  bullet  separates 
from  the  muzzle  is,  in  our  non-dimensional  time,  practically 
equal  to  10.  The  position  of  the  recompression  shock  is  shown, 
at  successive  instants  of  time,  in  Fig.  12.  Here  we  may  note  an 
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interesting  numerical  effect,  which  confirms  our  belief  that 
3hocks  should  be  fitted.  In  the  initial  phase  of  evolution  the 
shock,  as  computed,  clearly  shows  two  parts,  one  of  which  is  what 
[20]  calls  a  Mach  disc  and  the  other  is  the  beginning  of  a  plume 
shock;  both  are  in  the  right  position,  according  to  the  experi¬ 
mental  evidence.  Unfortunately,  the  way  the  imbedded  shock  is 
computed  does  not  allow  a  single  shock  with  a  sharp  corner  to 
subsist  indefinitely;  the  corner  itself  is  dragged  backwards  by 
the  'plume'  branch  which,  in  turn,  develops  oscillations,  sending 
numerical  errors  to  other  regions  of  the  flow.  Since  a  perfect 
determination  of  the  plume  shock  was  not  considered  of  primary 
importance  in  this  phase  of  the  research,  that  branch  was  artifi¬ 
cially  deleted  and  not  allowed  to  form  again.  In  other  words, 
our  current  results  reflect  an  attempt  to  'capture'  the  plume 
shock;  the  effect  is  clearly  shown  in  Fig.  12;  the  edge  of  the 
plume  is  pushed  outwards  and  the  Mach  disc  (which  itself,  for  be¬ 
ing  fitted,  is  properly  located)  is  much  wider  than  in  the  exper- 
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iments. 


To  complete  the  description  of  the  flow,  we  present  a 
streamline  picture,  a  plot  of  constant  Mach  lines,  a  plot  of 
isentropics  and  a  plot  of  velocity  vectors  in  Figs.  13  and  14. 
All  these  figures  refer  to  the  instant  of  emersion  of  the  bullet. 


15 .  Second  phase  of  precur sor  shock  evolution 


The  last  output  of  the  run  presented  in  the  preceding 
section  shows  the  perturbed  flow  field  around  the  muzzle  a3  the 
bullet  is  about  to  emerge  from  the  barrel.  To  continue  the 
analysis  until  the  bullet  is  totally  emerged,  the  computational 
program  must  be  modified.  Although  only  few  statements  have  to 
be  added  or  changed  in  the  program  used  for  the  preceding 
analysis  in  order  to  produce  the  program  used  for  the  current 
one, the  underlying  algebraic  manipulations  are  much  more  compli¬ 
cated  because  the  X=0  line  is  now  the  image  of  the  bullet's  mov¬ 
ing  wall  (in  full  or  partially);  consequently,  a  time-dependent 
mapping  function  is  to  be  used. 

In  lieu  of  (22),  we  use  the  mapping  sequence: 

z  =  (r^*)  [  (z2-1/z2)/2  -  logz2  -  i"] 

z1  ♦  1/z1  =  2B(C2  +  1/t2) 

(;2-  1  )/(C2+  1  )  =  [  (C.|-  D/O^f 

C1  =  (C  -  d2/0/d  -  d2) 

where  3  and  d  are  functions  of  time.  This  means  that  the  time 
associated  with  the  coordinates,  x  and  y,  of  the  physical  plane, 
cannot  be  confused  with  the  time  associated  with  the  coordinates, 
o  and  9  ,  of  the  mapped  plane.  We  will,  thus,  denote  the  first 
by  t  and  the  second  by  The  value  of  d  is  determined  to  assure 
correspondence  between  C  =  i  and  z  =  3,  where  3  is  the  abscissa 
of  the  nose  of  the  bullet.  The  value  of  3  is  made  equal  to  1/2 
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when  the  cylindrical  part  of  the  bullet  begins  showing  out  of  the 
barrel,  and  kept  constant  thereafter.  During  the  emersion  of  the 
bullet  nose,  6  is  varied  between  0  and  1/2  using  some  arbitrary 
law,  for  example  linearly. 


The  new  mapping  is  provided  to  maintain  the  computational 
grid  as  close  to  orthogonal  as  possible,  despite  the  appearance 
of  a  new  boundary,  and  to  keep  the  line  X=0  as  close  as  possible 
to  a  P=constant  line.  By  no  means,  however,  can  such  a  simple 
mapping  as  (89)  make  the  X=0  line  to  coincide  exactly  with  a 
Psconstant  line.  It  is  necessary,  thus,  to  introduce  the  image 
of  the  boundary,  opposite  to  the  precursor  shock,  in  the  form 
p  =  b(8,T),  and  to  normalize  the  grid  with  the  following 
transformation,  which  replaces  (48): 


X  =  (P  -  b)/(c  -  b) 

*  _  9  . tanh[«(Y-1 ) ]  * 

2  ”  o;  tanha  *  2 

T  =  t 


(90) 


In  addition  to  the  definitions  of  z,  C,  g,  ♦,  t  and  2  given  by 
(20),  (21),  (23),  (26)  and  (24),  we  need  to  define  the  quanti¬ 

ties: 


f  *  iT?-i  *  f,  *  ‘  f2 

•  •  ■  *1  *  1 


From  (89)  we  obtain 


f  =  T 


2d  1 ,  J  V1  ,  V1 

?(;,-  T)d. - —  log  - - r  « 


1-d 


21  ;  t 


26 


t 


dt 
dt , 


(9D 

(92) 


(93) 


*  s  "  2dV 


1-d 


.  2  .2 
C  «-d 


]  ♦  2f 


2  2 
C  *d 


6  [1  *  S 


S'1 

log  - 7  ]5 


’1 


(94) 


The  relations  given  by  (28),  (29),  (30)  and  ( 3 i )  must  be  comple¬ 
mented  by  the  following  expressions: 
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Pt  *  Pf1  •  9t  =  f2  *  %  -  °  *  Ty  =  0  •  Tt  *  1  (95) 

xT  =  (Sf^f^P/G  ,  yT  =  -(2f1+ef2)P/G 

(96) 

tp  =  0  ,  tfl  =  0  ,  tT  =  1 

Gx  =  8(*1-41f1*»2f2)  .  *T  =  *24,Wi  <97) 

Between  p,  9,  T  and  X,  Y,  T,  the  following  relations  hold: 

°X  *  1/Xp  •  PY  =  -  V(XPV  •  PT  =  -  XT/XP 

9X  =  o  ,  ey  s  t/Ye  ,  0?  =  o  (98) 

Tx  =  0  ,  Ty  =  0  ,  xT  =  1 

xp  =  1/PX  ,  X9  =  -  py/(px9Y)  •  XT  =  ■  °T/Px 

Yp  =  0  ,  Ye  =  1/ey  ,  YT  =  0  (99) 

Tp  =  0  ,  T0  =  0  .  Tt  =  1 

Consequently, 

xx  =  e/GXp  ,  x y  =  -(ZX9/pXp+S)p/GY9  ,  xT  =  -ZXT/GXp-(Zf 1-Sf2)p/G 
yx  =  a/GXp  ,  yy  =  (e-3X9/PXp)P/GY9  ,  yT  =  -SXT/GXp-(2f ^f^P/G 
tx  =  0  ,  ty  =  0  ,  tT  =  1  (100) 

and 

Gx  =  (Z^^*2)G2/p  ,  Gy  =  (2«1-«^2)G2/p  ,  Gt  =  G*1 

1  c  "  (101) 
“x  s  (Z*2-2*1)G/p  ,  wy  =  (Z*1+S«2)G/p  •  =  *2 

From  (90),  it  follows  that 

Xp=1/(c-b) ,  X9xXp((X-1)be-Xc9J,  XT=Xp[(X-1 )bT-XcT]  (102) 
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9y  =  [D2  -  (9-ir/2)2]a/D 


(103) 


with 


D  =  (*  -  9  )/tanh  “ 
o 


(104) 


Certain  second  derivatives  are  also  needed,  to  evaluate  viscous 
terms : 


xpp=°-  Xpe=Xp(be~c9) 


xee=xp9x0/xp+xp  ^xe  (b0-c0)-*-(x-i  )b8e-Xc0e] 


(105) 


Ypp=0  •  Ype=°  •  Yee  =  (2a/D)  Ye  (9-*/2> 


Let  Q  be  a  point  on  the  physical  plane: 


Q  s  xl  +  y3 


(106) 


Because  of  (100),  it  turns  out  that 


G  Qx  s  (1/Xp)I  ,  GYgQy  =  -  (X0/Xp)i  ♦  pj 
G  qt  =  -(XT/Xp>Pf1)l  -  Pfp 
The  normal,  S  to  an  X  =  constant  line  is  defined  by 


(107) 


f»  .  qy  =  0 


(108) 


Therefore,  its  components  in  the  I -and  j-directions,  respective¬ 
ly,  are  • 


N1  =  1/v  ,  N2  =  (X0/PXp)N7 


(109) 


with 


V  =  [U(Xg/PXp)2]1/2 


(110) 
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and 


At  the  body,  o  =  b,  X0/PXp  =  -b0/b,  XT/Xp  =  -bT ,  bT  = 


bt  = 


bT  +  Pt 


b09t  = 


br  +  pfi  +  bef2 


(ill) 


Let 


F(x,y.t)  =  0  (112) 

be  the  equation  of  the  X  =  0  line  in  the  physical  plane.  Then 
Fy  =  0,  that  is 


F  (Zbfl-Sb)  +  F  (2bfl+Zb)  =  0  (113) 

X  o  y  H 

Therefore , 

b„/b  =  (2F  -ZF  )/(ZF  *2F  )  (114) 

9  x  y  x  y 

Similarly,  F^.  =  0,  that  is 

Fx[Zbx-b(ef1-«f2)]+Fy[2bT-b(2f1+Zf2)]+OFt  r  0  (115) 

Therefore , 


b[(€f1-2f2)Fx+(2f1+Zf2)Fy]-GFt 

bt  =  ZF  +2F 

x  y 


(116) 


In  the  present  run,  the  bullet  has  been  assumed  as  a 

cylinder,  with  an  elliptic  nose.  Let  x  be  the  abscissa  of  the 

o 

center  of  the  ellipse,  which  travels  along  the  centerline  at  a 
constant  3peed;  let  a  be  the  major  axis  of  the  ellipse.  There¬ 
fore,  for  the  nose. 


_  ,  ,2  2  2  2  . 
F  =  (x+x  )  +a  y  -a  =0 
o 


(117) 


fx  =  2(x+xq)  ,  f  =  2a  y  ,  Ffc  =  2(x+xo)xQt  (118) 


In  the  cylindrical  part 
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F  =  y+1  ,  F  *  0  ,  F  =  1  ,  F  =  0  (119) 
x  y  u 

The  calculation  is  then  carried  on  in  a  way  similar  to 
the  one  performed  in  the  absence  of  the  bullet.  Obviously,  the 
velocity  of  the  flow  on  the  wall  of  the  bullet  equals  xqI.  In 
addition,  as  a  first  step  towards  the  simulation  of  spillage 
between  bullet  and  barrel,  the  pressure  at  the  rim  of  the  bore  is 
allowed  to  increase  linearly  with  time,  beginning  when  the  nose 
of  the  bullet  is  almost  completely  out. 


A  further  remark  is  necessary  for  a  proper  evaluation  of 
the  shocks.  The  image  of  a  shock  in  the  mapped  plane  is 


In  the  (P,9)  plane,  we  see  variations  in  the  time  t ; 


•dC  dz  31ogC 

(c  . _ s  _ s  3 

s;t  '  dz  dt  +  s  3t 


(121) 


The  first  term  in  the  right  hand  side  is  produced  by  the  varia¬ 
tion  of  za  in  time,  the  second  term  by  the  changes  in  the  mapping 
parameters  (the  derivative  is  made  at  a  constant  zs).  The  3hock 
point,  however,  is  always  taken  on  a  8  =  constant  line.  There¬ 
fore  , 


<y,  •  <v. 


i6 

e  =  (o  )  c  /p 

s  r  s  s 


(122) 


that  i3. 


(p  )TC  /p  =  g  dz  /dt  +  C  f  (123) 

s  T  3  3  “  3  3 

With  Q  defined  by  (106), 


W  =  5?.  HI  (124) 

at 

(33)  can  be  used  to  replace  I  and  J  and  then  (123)  to  replace  x 
and  y  ;  after  simplifying,  one  obtains 

3v 

(P  )  =  GW/N.  +  c  (f.  +  f_N_/N.)  (125) 

3  T  1  I  d  d  1 
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which  replaces  (86). 

The  sequence  of  isobar  plots  of  Fig.  15  describes  the 
computed  flow  field  around  a  protruding  bullet  which  has  the 
speed  and  the  size,  if  not  exacted  the  shape,  of  the  bullet  in 
[20],  Contrary  to  all  the  preceding  plots,  the  ones  in  Fig.  15 
are  confined  to  a  relatively  Small  region  in  the  vicinity  of  the 
bullet,  since  the  rest  of  the  flow  is  uneventful.  The  mesh  is 
indeed  very  coarse;  nevertheless,  a  good  amount  of  information  is 
still  obtainable.  / 

/ 

It  should  be/noted  that  the  speed  of  the  flow  in  front  of 
the  bullet  is  h^igher  than  the  speed  of  the  bullet  itself;  in 
fact,  the  flow  i j/  accelerated  by  the  supersonic  expansion  occur¬ 
ring  around  thf'lip  of  the  muzzle.  Therefore,  the  advancing  bul¬ 
let  acts  on  the  surrounding  flow  as  a  receding  wall  acts  on  a  gas 
at  rest,  and  produces  a  further  expansion. 


16. /Conclusions  and  recommendations 


The  results  of  our  numerical  tests  show  that  the  basic 
features  of  our  technique,  viz.  grid  definition,  integration 
scheme  and  3hock  fitting,  do  indeed  provide  very  high  accuracy. 
The  minimal  amount  of  grid  points  used  also  assures  that  the 
present  calculations  are  much  faster  than  any  other,  provided 
that  the  stepsizes  are  comparable.  No  total  running  time  on  a 
large-scale  computer  can  be  mentioned  yet,  since  all  calculations 
have  been  performed  on  a  minicomputer  (PDP  11/40)  which  not  only 
is  very  slow,  by  current  standards,  but  also  has  storage  require¬ 
ments  that  in  turn  force  the  program  to  be  written  in  a  highly 
inefficient  way. 

Before  giving  the  analysis  a  full  certificate  of  relia¬ 
bility,  however,  further  research  is  needed.  The  shock-fitting 
section  of  the  program  must  be  recoded,  in  order  to  let  the  plume 
branch  of  the  shock  be  tracked;  this  may  imply  a  deep  revision  of 
some  of  the  discretized  equations  used  in  the  fitting  because  the 
plune  shock  stretches  obliquely  with  respect  to  the  computational 
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grid,  whereas  the  recompression  shock  in  the  vicinity  of  the 
centerline  tends  to  lie  parallel  to  orconstant  lines.  In  addi¬ 
tion,  some  numerical  difficulties  may  arise  at  the  Intersection 
of  the  two  branches  of  the  shock,  which  has  the  physical  nature 
of  a  triple  point. 

The  dependence  of  the  results  on  the  Reynolds  number 
should  also  be  investigated  thoroughly.  So  far,  we  have  concen¬ 
trated  grid  lines  near  the  rigid  wall  but  we  have  made  no  attempt 
to  compare  our  flows  to  boundary  value  results  in  that  neighbor¬ 
hood.  Moreover,  the  edge  of  the  plume  i3  a  region  of  high 
viscous  stresses,  which  our  mesh  is  totally  inadequate  to  calcu¬ 
late.  Three  different  approaches  should  be  considered,  to  im¬ 
prove  the  present  results.  The  first  consists  of  using  a 
stretching  in  the  9-direction  which  concentrates  grid  lines  where 
the  vortlcity  is  the  highest.  The  second  (and  more  difficult, 
but  perhaps  richer  in  consequences)  consists  of  a  fitting  of  the 
shear  layer,  the  middle  line  of  which  should  be  tracked  as  a  line 
floating  within  mesh  points. The  third  approach  would  pu3h  this 
concept  to  its  extreme  consequences,  totally  eliminating  the  cal¬ 
culation  of  viscous  terms  and  replacing  shear  layer  by  vortical 
discontinuities.  An  inviscid  calculation,  not  requiring  extreme¬ 
ly  fine  grids  and  the  processing  of  the  viscous  terms,  would 
yield  a  very  welcome  reduction  in  running  time. 

If  such  a  stage  could  be  reached,  we  would  consider  it 
proper  to  compare  results  with  those  of  a  (not  yet  existing) 
purely  inviscid  code,  where  the  mechanism  for  the  onset  of  a 
plume  configuration  is  provided  by  the  appearance  of  a  3hort 
shock,  normal  to  the  wall,  totally  independent  on  the  recompres¬ 
sion  shock,  as  described  in  the  present  paper. 
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